
;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   Interpolating to specified pressure levels

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
;  a = addfile("../wrfout_d01_2000-01-24_12:00:00.nc","r")
;  input_file = systemfunc("ls /public/users/wrf/wrfout/2012061900/wrfout_d01_*")

  a = addfile(input_file + ".nc", "r")

  wks_type = "png"
  wks_type@wkWidth = 800
  wks_type@wkHeight = 800

; Set some Basic Plot options
  res = True
  res@MainTitle                   = "WRF Pressure Level Contour"
  res@Footer = False
  res@NoHeaderFooter = True

  pltres = True
  mpres = True
  mpres@mpDataBaseVersion       = "MediumRes"
  mpres@mpDataSetName           = "Earth..4"
  mpres@mpGeophysicalLineColor      = "Black"
  mpres@mpNationalLineColor         = "Black"
  mpres@mpUSStateLineColor          = "Black"
  mpres@mpGridLineColor             = "Black"
  mpres@mpLimbLineColor             = "Black"
  mpres@mpPerimLineColor            = "Black"
  mpres@mpGeophysicalLineThicknessF = 2.0
  mpres@mpGridLineThicknessF        = 1.0
  mpres@mpLimbLineThicknessF        = 1.0
  mpres@mpNationalLineThicknessF    = 2.0
  mpres@mpUSStateLineThicknessF     = 1.0
  mpres@mpOutlineBoundarySets       = (/"National","Continental", "Geophysical"/)   ; turn on country boundaries
  mpres@mpOutlineSpecifiers = "China:states"
  mpres@mpGeophysicalLineColor      = "Black"       ; color of cont. outlines
  mpres@mpGeophysicalLineThicknessF = 2          ; thickness of outlines
  mpres@mpNationalLineThicknessF = 2          ; thickness of outlines


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  ; Reading times
  char_times = a->Times
  hours = WRF_Times2Udunits_c(char_times, 0)
  ntimes = dimsizes(hours)         ; number of times in the file
  wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)
  ; convert to Beijing Time
  hours = hours + 8
  times_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")
  timesTideStr_cst = ut_string(hours, "%Y%N%D%H%M%S")
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file

; The specific pressure levels that we want the data interpolated to.
;  pressure_levels = (/ 850., 700., 500., 300./)   ; pressure levels to plot
  pressure_levels = (/ 500./)   ; pressure levels to plot
  nlevels         = dimsizes(pressure_levels)     ; number of pressure levels

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  contour_dt = 1 ; time interval of contour
  dt = tointeger(contour_dt / wrf_dt)
  do it = 0,ntimes-1,dt             ; TIME LOOP

    output_file = output_prefix + timesTideStr_cst(it)
    wks = gsn_open_wks(wks_type, output_file)

    print("Working on time: " + times_cst(it) )
    res@TimeLabel = times_cst(it)   ; Set Valid time to use on plots

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        

    tc = wrf_user_getvar(a,"tc",it)        ; T in C
    u  = wrf_user_getvar(a,"ua",it)        ; u averaged to mass points
    v  = wrf_user_getvar(a,"va",it)        ; v averaged to mass points
    p  = wrf_user_getvar(a, "pressure",it) ; pressure is our vertical coordinate
    z  = wrf_user_getvar(a, "z",it)        ; grid point height
    rh = wrf_user_getvar(a,"rh",it)        ; relative humidity

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    do level = 0,nlevels-1                 ; LOOP OVER LEVELS

      pressure = pressure_levels(level)

      tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)
      z_plane  = wrf_user_intrp3d( z,p,"h",pressure,0.,False)
      wrf_smooth_2d(z_plane, 3)
      rh_plane = wrf_user_intrp3d(rh,p,"h",pressure,0.,False)
      u_plane  = wrf_user_intrp3d( u,p,"h",pressure,0.,False)
      v_plane  = wrf_user_intrp3d( v,p,"h",pressure,0.,False)

      spd     = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec
      spd@description = "Wind Speed"
      spd@units = "m/s"
      u_plane = u_plane*1.94386     ; kts
      v_plane = v_plane*1.94386     ; kts
      u_plane@units = "kts"
      v_plane@units = "kts"


      ; Plotting options for T                
;        opts = res                          
;        opts@cnLineColor = "Red"
;        opts@ContourParameters = (/ 1.0 /)
;        opts@cnInfoLabelOrthogonalPosF = 0.07  ; offset second label information
;        opts@gsnContourLineThicknessesScale = 1.5
;        contour_tc = wrf_contour(a,wks,tc_plane,opts)
;        delete(opts)


      ; Plotting options for RH                
;        opts = res                          
;        opts@cnFillOn = True  
;        opts@pmLabelBarOrthogonalPosF = -0.1
;        opts@ContourParameters = (/ 10., 90., 10./)
;        opts@cnFillColors = (/"White","White","White", \
;                              "White","Chartreuse","Green",\
;                              "Green3","Green4", \
;                              "ForestGreen","PaleGreen4"/)
;        contour_rh = wrf_contour(a,wks,rh_plane,opts)
;        delete(opts)


      ; Plotting options for Wind Speed                
;        opts = res                          
;        opts@cnLineColor = "MediumSeaGreen"
;        opts@ContourParameters = (/ 10. /)
;        opts@cnInfoLabelOrthogonalPosF = 0.07  ; offset second label information
;        opts@gsnContourLineThicknessesScale = 2.0
;        contour_spd = wrf_contour(a,wks,spd,opts)
;        delete(opts)


      ; Plotting options for Wind Vectors                 
        opts = res          
        opts@FieldTitle = "Wind"   ; overwrite Field Title
        opts@NumVectors = 47       ; wind barb density
        vector = wrf_vector(a,wks,u_plane,v_plane,opts)
        delete(opts)


      ; Plotting options for Geopotential Heigh
        opts_z = res                          
        opts_z@cnLineColor = "Blue"
        opts_z@gsnContourLineThicknessesScale = 2.0


      ; MAKE PLOTS                                       


        if ( pressure .eq. 500 ) then   ; plot temp, height, wind barbs
          opts_z@ContourParameters = (/ 50.0 /)
          contour_height = wrf_contour(a,wks, z_plane,opts_z)
          plot = wrf_map_overlays(a,wks,(/contour_height, \
                                    vector/),pltres,mpres)
        end if

        delete(opts_z)
        delete(wks)

    end do      ; END OF LEVEL LOOP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    

  end do        ; END OF TIME LOOP

end
